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We construct a 3-color Potts-like model for the graphene zigzag edge reconstructed with Stone- 
Wales defects, in order to study its thermal equilibrium properties. We consider the cases in which 
the edge's dangling carbon bonds are both non-passivated or fully passivated with hydrogen. We 
show that in both cases there is always a finite concentration of defects at any finite temperature 
and moreover, that such concentration is exponentially dependent on the effective parameters that 
describe the model, which were estimated using DFT results both from the literature, as well as our 
own. Such equilibrium mechanisms place a lower bound on the concentration of defects in zigzag 
edges, since the formation of such defects is due to non-equilibrium kinetic mechanisms. 

PACS numbers: 65.80.Ck, 72.80.Vp, 05.50.+q 



I. INTRODUCTION 

The mechanical exfoliation of graphene,^ a one-atom 
thick sheet of carbon atoms arranged in a honeycomb 
structure, followed by subsequent experiments E1S1 nas 
given rise to a profusion of studies, both experimental 
and theoretical (see^l for a list of references) regarding 
the physical properties of this remarkable material. Such 
an interest has not subsided to the present date, quite on 
the contrary. 

A simple nearest-neighbor tight-binding approxima- 
tion of the electronic Hamiltonian in graphene reveals 
that the honeycomb lattice structure leads to a disper- 
sion relation that is linear around two specific points of 
the Brillouin zone, the Dirac points. Since the Fermi 
level of pristine graphene lies at these points, its quasi- 
particles behave, in a continuum approximation, as mass- 
less relativistic fermions with a speed of light equal to the 
Fermi-velocity (vf ~ lC^rns -1 )! 5 -^ 

Typically, graphene exhibits high crystal quality, large 
transparency, being highly conductive and very strong 
yet flexible. 5 6 All these characteristics place graphene as 
a good candidate to be used in a variety of technological 
applications, such as in solar cell technology,^ in liq- 
uid crystal devices,^ in single molecule sensors,^ in the 
fabrication of nano-sized prototype transistors™ among 
many others. Understanding the transport properties of 
graphene is thus an essential research program towards 
its application to future nanoscopic devices. 

Several of these sought nanoscopic devices will cer- 
tainly be based on the use of graphene ribbons and 
graphene quantum dots. In particular, graphene ribbons 
are usually classified as zigzag or armchair, depending 



on their edge configuration (see Fig. |T|). It is already 
well established that the electronic properties of these 
nanostructures are strongly affected by their edge con- 
figuration. A nearest-neighbor tight-binding approach 
leads to the conclusion that zigzag ribbons are metallic 
irregardless of their width, while armchair ribbons can 
be either semiconducting or metallic, depending on their 
width. In addition, zigzag ribbons present edge local- 
ized states around the Fermi energyPSHn] However, ab- 
initio calculations predict that graphene ribbons are al- 
ways semiconducting.^ Experimental results show that 
the ribbons' energy gaps increase with decreasing ribbon 
widthP 
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Figure 1. (a) Scheme of zigzag ribbon, (b) Scheme of an 
armchair ribbon. 

It has been showrJ 2 ^!! that edge disorder modi- 
fies substantially the electronic properties of nanostruc- 
tures based on graphene and is responsible for most of 
the transport properties in these systems. The wave- 
functions associated with zigzag edges and vacancies^! 
decay very slowly with the distance from the edge be- 
cause of the absence of a gap in the graphene spectrum. 
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This slow decay leads to strong quantum interference ef- 
fects that are responsible for destructive quantum inter- 
ference, Anderson localization, and Coulomb blockade 
effects in the electrical conductance of graphene-based 
devices While most of edge disorder in graphene nanos- 
tructures is produced by the method of cutting graphene 
by hot plasmapHSD it a i so shows conspicuously in chem- 
ically driven methods which can be classified as quasi- 
cquilibrium. 27 It thus follows that equilibrium mecha- 
nisms, such as the ones described in this work, only place 
a lower bound on the amount of disorder that can exist 
in these nanostructures and hence an upper limit in the 
value of the conductance that can be obtained in these 
devices. 

It is thus of extreme importance for graphene electron- 
ics development to understand the effect of edge disor- 
der in the transport properties of graphene ribbons and 
quantum dots. 

One example of defects appearing at graphene edges 
are sets of five and seven sided rings of carbons, usually 
named Stone- Wales (SW) defects 28 . These defects have 
been observed and found to be meta-stable in the bulk of 
graphene sheets. 21 ' Moreover, and regarding graphene's 
high temperature behavior, it was found that the forma- 
tion of Stone- Wales defects is the first step in the process 
of graphene's melting (T me i ting s» 4900K). 30 Besides, it 
has been shown through ab-initio calculations that when- 
ever SW defects are present in non-passivated graphene 
nanoribbons, the ribbons energy decreases as the defect 
approaches the edge of the ribbonP^In addition, further 
ab-initio calculations have also shown that the forma- 
tion of SW defects at the edges of both armchair and 
zigzag nanoribbons, stabilize them, both energetically 
and mechanicallyp2H2lln particular, in the absence of hy- 
drogen passivation, the zigzag edge is meta-stable under 
total reconstruction with SW defects! 23 ' Moreover, such 
total reconstruction of the zigzag edge gives rise to the 
appearance of a new kind of edge statePSESl However, if 
the zigzag edges are hydrogen passivated, it is the perfect 
zigzag edge that has the lower energy. The reconstruc- 
tion of the zigzag edge by SW defects acts as a mecha- 
nism that self-passivates the edgeP^l Density functional 
theory and molecular dynamics calculations corroborate 
these results, pointing to an energy barrier associated 
with the edge reconstruction of about 0.4-0.9eV per edge 
unit C eU.IEEZIBg 

These types of zigzag edge reconstructions (see Fig. 
|2|, are claimed to be stable only at very low hydrogen 
pressure (well below ambient conditions) and very low 
temperatures! 4 ^ However, reconstructions of the zigzag 
(as well as armchair) edges have been recently observed 
with high-resolution TEM,EilE2l albeit under rather ex- 
treme conditions, namely, the graphene flake is bom- 
barded with high-energy electrons (80keV) that remove 
C atoms from the sheet. The recent work of Suenaga 
et aZ.P on single-atom spectroscopy using low- voltage 
STEM, provides a non-destructive method of identifying 
the edge configuration of graphene ribbons, as does the 



work of Warner et al. on the observation of real-time 
dynamics of dislocations using high-resolution TEM. 45 
Moreover, refinements in other techniques, such as Ra- 
man spectra of the edges STM images of the edgesp^l 
or coherent electron focusing^! may help in the identifi- 
cation of edge reconstructions. 
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Figure 2. (a) Scheme of clean zigzag edge, usually named zz. 
(b) Scheme of a totally reconstructed (with SW defects) zigzag 
ribbon, usually named zz(57). In (c) and (d), we present 
schemes of a partially reconstructed zigzag ribbon, respec- 
tively, zz(576) and zz(5766). 

In this work, we construct a one-dimensional 3-color 
Potts-like model so as to describe the reconstruction of 
the zigzag edge with SW defects, due to thermal fluc- 
tuations. From such a model, we will extract thermo- 
dynamic properties of the zigzag edge, both in the ab- 
sence and in the presence of hydrogen passivation, such 
as the number of pentagons and heptagons in the edge, 
the density of links between pentagons or heptagons and 
hexagons, the density of domains formed between groups 
of pentagon-heptagon pairs, groups of hexagons, the av- 
erage size of these domains and the domain size distribu- 
tion of domains of each kind at the edge. These quantities 
are then used to characterize the degree of disorder of the 
zigzag edge due to thermal fluctuations. 

The simple model presented in this work indicates that 
if the nanoribbon edge behaves as a one-dimensional sys- 
tem with short-range interactions, it will always present 
a finite concentration of defects at any finite temperature 
and moreover that such a concentration is exponentially 
dependent on the values of the effective parameters that 
enter the model, which we are able to extract from DFT 
calculations. However, the lower the temperature, the 
scarcer and the smaller the defective domains will be. 
Here, we name as defective domains those that have a 
different spin-order from that of the ground state. 

The structure of this paper is as follows: Section [TT] will 
be devoted to the presentation of the 3-color Potts-like 
model that describes the thermodynamics of the edge 
with SW defects and of the results extracted from it. 
We will first present a descriptive outline of the model. 
In sub-Section |II A| we will show how the exchange pa- 
rameters that enter the model can be extracted from the 
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ab-initio results. In sub-Section |II B[ we will compute the 
thermodynamic quantities characterizing the edge, using 
a transfer matrix formulation of the Potts- like model and 
we will analyze their dependence both on the tempera- 
ture and on the exchange parameters. Finally in section 
|III| we will present our conclusions. We leave to the 
appendices the actual numerical computation of the ex- 
change parameters from particular ab-initio results, the 
calculation of correlation functions in the 3-color Potts- 
like model, using the transfer matrix formalism, and the 
explicit computation of the domain size distribution of 
domains of polarized and unpolarized spins. 



II. POTTS-LIKE MODEL OF THE ZIGZAG 
EDGE 

We will study graphene zigzag edges with SW defects, 
employing a one-dimensional 3-color Potts-like model, 
where each color is assigned to a different polygon of 
the edge (hexagons, heptagons and pentagons). We la- 
bel each edge unit cell, i. e. each polygon at the edge, 
by the integer variable i = 0, 1, ■ • • , 2N, with the state 
of such a cell being described by the ternary variable 
<7j = 0, ±1, according to whether the polygon forming 
that cell in the reconstructed edge is an hexagon, hep- 
tagon or pentagon. We consider a nearest-neighbor cou- 
pling between adjacent cells only (the validity of this as- 
sumption will be justified in Appendix [A]) , which leaves 
us with 9 possible values for the couplings J CTi(Ti+1 , de- 
pending on the neighboring states. We take as refer- 
ence state with zero energy the perfect zigzag edge, thus 
Jqo = 0. Taking into account the experimental observa- 
tions, we will exclude from the model states where two 
pentagons or heptagons sit at neighboring sites, i.e. pair- 
ings of heptagons or pentagons are forbidden and one 
has J ++ = J = oo. Invariance under inversion im- 
plies that the order in a pentagon-heptagon, pentagon- 
hexagon or heptagon-hexagon pair is irrelevant, and thus 

J |_ = J_| , Jo+ = J+o and Jo- = J-q. Moreover, since 

heptagons and pentagons are created in pairs through 
the transference of C atoms between neighboring sites, 
we will assume that the probability of creation of a pen- 
tagon or an heptagon is the same, which implies that 
Jo- = Jo+- Hence, the 9 initial possible values of the cou- 
plings are reduced to two free parameters: Jo+ = 7 > 0, 
which reflects the fact that the formation of defects costs 

energy and J_| = S, which may be negative or positive 

depending on whether the totally reconstructed edge has 
lower or higher energy than the pristine zigzag edge (i.e., 
depending on the state of passivation of the edges, as 
discussed above). Finally, since C atoms are conserved, 
one should, strictly speaking, consider a model with as 
many heptagons as pentagons, i.e. one should work in 
a subspace of the state-space having the overall mag- 
netization M — Y^i=o °~i = 0- Such a constraint can 
be written in terms of an imaginary applied magnetic 
field over which one has to integrate, once the eigenval- 



ues of the transfer matrix of the Potts-like model have 
been computed. In such a case, the eigenvalues can no 
longer be simply determined. We will therefore relax 
this constraint and we will only implement it on aver- 
age, as ( M ) =0 in Id. Note, moreover^that although 
some of the edge observation techniques"^ are highly 
energetic and cause the ejection of C atoms from the 
edges, the system cannot be considered to be in thermo- 
dynamic equilibrium when such ejection occurs and the 
model introduced below is therefore not applicable.^ It 
may however be applicable after a characteristic relax- 
ation time, such that the thermodynamics of the edge 
would be described in terms of an effective temperature, 
dependent on the energy deposited by the electron beam 
and the heat conduction process in graphene. One would 
expect the number of (remaining) C atoms in the edge 
to be conserved in this late-time regime. In figure [3] we 
present a cartoon of three possible configurations of the 
edges and how they translate into configurations of the 
3-color Potts-like model. 



A. The exchange parameters from ab-initio results 

In what follows we will show how to compute the ex- 
change parameters of the Potts-like model, from ab-initio 
results of zigzag ribbons with reconstructed edges. 

The energy per edge atom of periodic edge configura- 
tions such as zz, zz(57), zz(576) and zz(576") (where n 
stands for the number of hexagons in the periodic edge 
configuration), can be readily computed either using the 
3-color Potts-like model proposed above (for particular 
values of the parameters 7 and S), or using the ab-initio 
results for the edge energies computed from density func- 
tional theory. From a least squares method, we can read- 
ily compute the exchange parameters, 7 and 5, of the 
Potts-like model, in such a way that the latter describes, 
to a good degree of accuracy, the ab-initio results. 

The edge energies (per unit cell of the perfect zigzag 
edge) of different reconstructed edges, for example, 
zz(57), zz(576), zz(5766), etc., are given, in the scope 
of the previously introduced Potts-like model, by 

E(zz{57)) = J+-, (la) 

E(zz(576 n )) = J +-+ 2J +" + ("- 1 ) J »" ; (lb) 

E(zz) = J 00 , (lc) 

where n stands for the number of edge hexagons present 
in a unit cell. Expressing these energies relative to 
the clean edge energy AE(zz(576 n )) = E(zz(576 n )) - 
E[zz), with the latter set to zero (i.e. J o = 0, as above), 
we obtain, 

AE(zz(57);S, 7) = 8, (2a) 

A£(zz(576"U 7 ) = ^J, (2b) 
n + 2 

where we made the substitutions J_| = S and J + o = 7- 
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Given the extreme sensitivity of the thermodynamic 
quantities on the values of exchange parameters 7 and 5, 
an example of the actual computation of these parame- 
ters is left to Appendix |A"| 

B. Edge thermodynamics 
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Figure 3. Scheme of the Potts-like model (3-color) for the 
zigzag edge with SW defects. In (a), a clean zigzag edge, zz, 
is shown. In (fa), a zigzag edge which is totally reconstructed 
with SW defects, zz(57), is presented. In (c), a zigzag edge 
with an arbitrary reconstruction is shown. 



We now consider the energy e n of the edge zz(576") 
referred to the pristine zigzag edge, as obtained from ab- 
initio calculations (see Fig. [9]). The exchange parame- 
ters 7 and 5 can be obtained from a minimization of the 
sum of the squared differences between A_E(zz(576 n )) , 
as given by Eqs. Q, and e„, 

n 2 



S(S,7) = J2 [A£(zz(576 n );J )7 ) 



(3) 



The uncertainty on the computed exchange parameters, 
is given by 



71=0 



/ dz 
\ de„ 



(4) 



where z stands for 7 or 6, whose expression as a func- 
tion of e n is computed from minimization of S(5, 7), and 
where o~ n are the uncertainties in the ab-initio energies 



In the previous paragraphs, we have shown how to 
map the different configurations of a reconstructed edge 
of a graphene zigzag ribbon to those of a 3-color nearest- 
neighbor Potts-like model. We have also shown how the 
two parameters that characterize the coupling between 
the neighboring spin variables can be extracted from ab- 
initio results. We now wish to use such a model to com- 
pute useful quantities relating to the thermodynamics of 
the edges. As is usual in one-dimensional models with 
nearest neighbor interactions, the speediest way to com- 
pute thermodynamic properties, including spin-spin cor- 
relation functions, is to express these quantities in terms 
of a transfer matrix!^ In the case of the model presented 
above, the transfer matrix reads: 



(5) 



where f3 — 1/fcgT and 7 and S are the parameters 
introduced in subsection |II A| The model, as defined 
by Eq. ([5|, represents a limiting case of the Blume- 
Emery-Grifhths model in one dimensionP2HHI One eigen- 
value of the transfer matrix, Ao = — e~@ s , can be read- 
ily identified, after which the other two are also easily 
computed from the quadratic equation that is obtained 
from the application of, e.g., Ruffini's rule to the cu- 
bic secular equation. These two eigenvalues are given 




by A± 
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At all 



-PS ± - e -psy + 8 e -2pj 

temperatures above zero, A+ is the largest eigenvalue. At 
T = and if S > 0, this is also the case, however if 5 < 
the largest eigenvalue may be doubly or thrice degener- 
ate, which reflects the degeneracy of the ground-state of 
the system (see Appendix [B| . 

The free energy of the system is given by F = 
—kgT ln(Tr T 2N ) from which we have that in the thermo- 
dynamic limit N — > 00 the free energy per site is simply 
proportional to the logarithm of the largest eigenvalue, 



-k B T\n{ - 



(6) 

We are primarily interested in the disorder caused ei- 
ther to a clean zigzag edge or to a totally reconstructed 
zigzag edge (zz(57) edge) through the effect of temper- 
ature, which leads these configurations [as depicted in 
Fig. |3ja) and (&)] to evolve into[3jc). A measure of such 
disorder can be obtained by counting the number of do- 
mains of polarized spins (heptagons and pentagons), that 
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Table I. Contribution of a H — domain for the different ob- 
servables. 



exist between sites with 0-spin (hexagons) , or vice- versa. 
As an example, in Fig. |3jc) one can count two domains 
of polarized spins. In order to be able to count them, 
consider the contribution of a domain, both to X^ ^ 
which measures the number of heptagons or pentagons 
in the system, and to <5 CTi(Ti+1 ._i, which measures the 
number of heptagon-pentagon links (see Table [T]). bmce 
each domain contributes exactly 1 to the difference be- 
tween these two quantities, one sees that the number of 
domains is given by the difference of these two opera- 
tors. Note however, that whenever the spin chain has no 
0-spins, the difference between these two operators gives 
0. As a consequence, the exact expression of the average 
domain density of ±-spin domains, (nd±) = (Nd±)/2N, 
is given by 



(n d± ) 



1 

2N 



J aiUi +1 , 



n« 



(7) 

When using periodic boundary conditions (PBCs), if 
the spin chain is not uniformly polarized (either all the 
sites with spin 0, or all the sites with polarized spin + 
or — ), the number of domains of ±-spins is always equal 
to the number of domains of 0-spins. In such a case, 
we have N d ± = N d o = N d . As a consequence, we can 
express (n d o) in terms of (n d ±), just by considering the 
following sum over all spin configurations, 
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.(8) 



Note that in the thermodynamic limit, 2N — > oo, the 
thermal averages of the products can be neglected, re- 
sulting in (n d0 ) « (n d ±). 

One can separately compute the correlation functions 
(of) and (S aiai+lt _i), as is done in Appendix [b| How- 
ever, it is simpler to consider instead generating fields 
in the partition sum that are coupled to ^ i of and to 
Y]j <5 CTi CTi + i ,—i ■ One then concludes, using the transfer 
matrix formalism, that the density of polarized sites, 

(n P oi) = %p = jn Si(°f )) is given, in the thermo- 
dynamic limit, by 
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85 



(9) 



The density of unpolarized sites is simply obtained from 
(n unp ) = 1 — (n po i). Moreover, the density of links 

between polarized sites, defined as, (n_| ) = 2iv = 

2^ 5tTi<Ti+i,-i )j is given, in the thermodynamic 
limit, by 

= 86 ■ 



(10) 



Similarly, the density of links between unpolarized sites, 

is simply obtained from (noo) = 1 — (n-\ ) — ( n ±o)- where 

n±o stands for the density of links between polarized and 
unpolarized sites. Note that in the thermodynamic limit, 
(?i ±0 ) s» 2(rido} ~ 2(n,d±)- As a consequence, in this 
limit, (n d ±) can be readily computed from Eq. (|9| and 
Eq. ( [To] ) . We can then write (i2d±) as 

2d 1 ' 

Substituting Eq. Q in Eq. ([9]), we obtain for (n po i) 

M = {1 + e - P 8 + e)g ( 12 ) 

where = >/(l - e~f 3& ) 2 + Se" 2 ^. A plot of this quan- 
tity as a function of T/j, for some values of the ratio 5/j, 
is shown in Fig. [4] 
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Figure 4. Plot of the density of polarized spins, (n po i), either 
for a non-passivated edge (in dark blue, violet and light blue) 
and for a hydrogen-passivated one (in red, orange and pink) 
as a function of T/j, for three different values of the ratio 
<5/7. The full (dark blue and red) lines stand for |<5/7| = 0.1; 
The dashed (violet and orange) lines stand for \5/^y\ — 1.; The 
dash-dotted (light blue and pink) lines stand for |<5/7| = 10.. 
Note that the S/y ratio is negative for the non-passivated case 
(because S < 0) and positive for the passivated case (because 
then (J > 0). The dashed green flat line represents the infinite 
temperature limit for both cases, which is 1/2. Since we are 
using PBC, {n unp } = 1 - (n po i). 



Substituting Eq. Q in Eq. ([ToJ, we obtain for (n 

. . e-"*(-l + e-^ + 0) 
(«+-) : 
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A plot of (13 1 as a function of T/j, for some values of 
the ratio Spy, is shown in Fig. [5] 
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Figure 5. Plot of the density of links between polarized spins, 
(n+_), either for a non-passivated edge (in dark blue, vio- 
let and light blue) and for a hydrogen-passivated one (in red, 
orange and pink) as a function of T/7, for three different 
values of the ratio (5/7. The full (dark blue and red) lines 
stand for | <5/'y j = 0.1; The dashed (violet and orange) lines 
stand for |<5/7| = 1-j The dash-dotted (light blue and pink) 
lines stand for | (5/-y | = 10.. Note that the 8/j ratio is nega- 
tive for the non-passivated case (because 8 < 0) and positive 
for the passivated case (because then S > 0). The dashed 
green flat line represents the infinite temperature limit for 
both cases, which is 1/(2 + 2^2). Since we are using PBC, 
(n o) = 1 - (n + _) - (n ±0 ). 



Finally, substituting Eq. 

(n<i±) 

(nd±) 



6| in Eq. (Ill, we obtain for 
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which gives the density of domains as a function of the 
temperature and of the coupling parameters. A plot of 
(nd±) as a function of T/7, is shown for some values of 
the ratio S/j in Fig. [6] 

In both the non-passivated case and the hydrogen- 
passivated one, the density of domains of polarized sites, 
rid, is very small at low temperatures (Fig. |6j. However, 
from Figs. [4] and [5] we conclude that these two situa- 
tions are substantially different. In the former, at low 
temperature, we have a small number of very large po- 
larized domains, with very few 0-spins between them. In 
contrast, in the latter case, at low temperature, we have 
a low number of very small polarized domains, with large 
domains of 0-spins between them. 

The domain size distribution (DSD) of domains with 
polarized spins, V±(L) — (N d t L /Nd±) , where L is the 
length of the domain and is the number of H — do- 
mains with size equal to L, can be computed exactly us- 
ing the transfer matrix formalism for the 3-states Potts- 
like model developed above. The detailed (and rather 
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Figure 6. Plot of the density of domains (of polarized spins), 
(n<i±), for a non-passivated edge (in dark blue, violet and 
light blue) and an hydrogen-passivated one (in red, orange 
and pink) as a function of T/7, for three different values of 
the ratio 5/7. The full (dark blue and red) lines stand for 
(5/7 = 0.1; The dashed (violet and orange) lines stand for 
<5/7 1 = 1.; The dash-dotted (light blue and pink) lines stand 
for |<5/7| = 10.. Note that the 5/7 ratio is negative for the non- 
passivated case (because 8 < 0) and positive for the passivated 
case (because 5 > in such a case) . The dashed green flat line 
represents the infinite temperature limit for both cases, which 
is l/(2+v / 2). Since we are using PBCs, in the thermodynamic 
limit, (n d o) ~ (w d ±). 



lengthy) calculation is presented in Appendix [C] We ob- 
tain, in the thermodynamic limit, the result 



P±(L) 



A + e 



p6 



(Ah 



B B8\ 



(15) 



Likewise, we have also computed the DSD of domains 
with unpolarized spins, T'o(L) = (N^/Ndo) (where AT° 
is the number of domains with size equal to L), see 
again Appendix [C] The result that we have obtained is 
given, in the thermodynamic limit, by 



A, 



1 



In Eqs. (151 and (16 1, A + = (1 



(16) 



0)/2. Equations 



( 15 I and ( 16 ) and their derivation are the main result of 
this work. Do note that L is geometrically distributed 
in both cases. In Fig. [7] we plot the DSD as a function 
of L (with logarithmic scale in the y-axis), for different 
values of T/7 and of the ratio We plot the DSD 

of domains of unpolarized spins when the edge is non- 
passivated [panel (a)] and the DSD of domains of polar- 
ized spins when the edge is hydrogen-passivated [panel 
(b)[- 

The characteristic functions of these two distributions 
can be readily computed from (15) and (16). We obtain 
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Figure 7. Plot of the DSD, P(L), for several temperatures and 
several values of the ratio 8/"/. The DSD is plotted with loga- 
rithmic scale in the y-axis. (a) DSD of domains of unpolarized 
spins in a non-passivated edge (7 > and S < 0). (b) DSD 
of domains of polarized spins in a hydrogen-passivated edge 
(7 > and 5 > 0). The red, blue, and green curves stand, 
respectively, for T/7 = 10000 K/eV, T/7 = 1000 K/eV and 
T/7 = 150 K/eV. The full and dashed lines, stand, respec- 
tively, for \S/j\ = 1.0 and \S/j\ = 0.1. 



in the case of domains with polarized spins, 

X+eP* - 1 



V±(w) 



\, e l3S e -i 



1 ' 



while for the case of unpolarized domains, we have 



Aj 



1 



(17) 



(18) 



The characteristic functions as given by Eqs. (17 1 and 



( 18 1 are analytic at w — 0, a result that implies that both 
distributions cannot have fat tails, since all the moments 
of these distributions exist. 

The first moment of these distributions, gives us the 
average size of, respectively, the domains of polarized and 
unpolarized spins. The explicit expression for the average 
size of the domains of polarized spins reads [see Appendix 
Eq. flC33l )], 



L+ = 



X+el 1 



\+e? s - 1 ' 



(19) 



while the average size of the domains of unpolarized spins 
reads 



L 



A+ 
A + - 1 



(20) 



It is interesting to compare the results given in ( 19 1 and 



( 20 1 with the results obtained from a different (and rather 



L± = (n po i) I (n d± ) and L 
for L±, the result 



(n unp ) / (n d0 ) . We obtain 



L+ = 



4e -2ff7 + e -^(_ 1 + e -^- + | 
4e-2/3 7 



(21) 



Moreover, in the thermodynamic limit, Lq can be writ- 
ten in terms of L± as L = (n unp ) / (ndo) = (l — 
{n po i}) I (n>do) ~ l/(n,do) — L±. If we now substitute A+ 



by its definition in equations (19 I and (20 I, we can show 
that L± = L± and Lq = Lo, i.e. the two definitions 
yield identical results. This equality suggests that the 
statistical variables Nd L /Nd and Nd are independent in 
the thermodynamic limit, i.e. that the fraction of do- 
mains with size L is independent of the number of the 
said domains in the limit of an infinite system. However, 
we were as yet unable to prove that such independence 
holds in the whole temperature range. 

In Fig. [5J we plot the average domain size of the minor 
domains in each of the two cases (Lo for non-passivated 
edges and L± for hydrogen-passivated edges - see the 
previous paragraph) as a function of T/8, for some values 
of the ratio 'y/S. 



Infinite temperature limit 
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Figure 8. Plot of the average domain size of the minority do- 
mains as a function of T/7. In dark blue, violet and light blue, 
the average domain size of unpolarized domain sites (unpassi- 
vated edge), Lo. In red, orange and pink, the average domain 
size of polarized domain sites (hydrogen-passivated edge) , L± . 
The different curves of the same kind, stand for three differ- 
ent values of the ratio 5/7. The full (dark blue and red) lines 
stand for |<5/7| = 0.1; The dashed (violet and orange) lines 
stand for | <5 / --yj = 1.0; The dotted (light blue and pink) lines 
stand for |<5/7j = 10.0. Note that the 5/7 ratio is negative 
for the non-passivated case (since S < 0) and positive for the 
passivated case (since S > 0). The dashed green flat line rep- 
resents the infinite temperature limit for both cases, which is 
1 + 1/V2. 



natural) definition of the average domain size, namely 



From Figs. |4][8] we confirm that irrespective of the 
value of the dimensionless temperature T/7, the sys- 
tem always presents a finite concentration of defects at 
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any finite temperature, as is to be expected for a one- 
dimensional system with short-range interactions. How- 
ever, since our model does not possess the full Z(S) sym- 
metry characteristic of a true Potts-like model (in which 
case the Peierls argument does apply), the results ob- 
tained are qualitatively different from those obtained for 
an Ising chain, where the formation of domains of macro- 
scopic size fully destroys order at any T — > (in our 
case, the energy of formation of a domain does depend 
on the domain's size). In contrast with what happens in 
the one-dimensional Ising model, here the disorder tends 
to zero with decreasing temperature, being exactly zero 
only at T = 0. Depending on the exchange parameters 
of the system, the concentration of minority domains at 
room temperature (and hence the degree of disorder of 
the edge at such temperature) may or may not present a 
large value. In addition, the smaller the ratio 7/(5 is, the 
less stable the edge is to the effect of thermal disorder. 
Furthermore, as expected, the average mean size of the 
minority domains increases with temperature. Finally, 
the larger the ratio 7/ J is, the larger is the size of the 
minority domains. This is to be expected, so as to mini- 
mize the number of domain walls (links 0+ and 0-, whose 
exchange parameter is 7) relatively to the number of +- 
links (whose exchange parameter is 8). 

In Appendix [AJ based on ab-initio calculations that 
we have both performed ourselves (case C\) and that 
we have obtained from the existing literature (case C2), 
we compute specific values of the exchange parameters 
for the introduced model, (<$Cd7Ci) an d (<5c 2 7 7c 2 )- Us- 
ing these particular values of the exchange parameters, 
in Appendix [X] we present plots of the thermodynamic 
quantities introduced in Eqs. (12l-(21l, as functions of 



the absolute temperature, rather than presenting them as 
functions of the reduced temperature T/7 and the ratio 
5/7. From these results, one can conclude that, in both 
particular cases C\ and C2 , the exchange parameters cal- 
culated from the ab initio are such that the ground state 
edge configuration is robust with respect to the effect of 
thermal disorder. 

Nevertheless, it should also be stated that due to the 
sensitivity of the system's thermodynamic behavior on 
the precise numerical value of the exchange parameters 
(see Appendix [Aj, this conclusion may well be challenged 
in the future, in case more detailed ab-initio calculations 
yield different numerical values for the exchange param- 
eters. We should emphasize in this regard the following 
facts: the ab-initio calculations we have performed, were 
done using narrow ribbons, where an interaction between 
the two edges of the ribbon can be observed; moreover, 
our ab-initio calculations did not take into account either 
the spin-polarization of electrons on the edge, or the re- 
laxation of the atoms along the transverse direction of the 
ribbon and that such complications need to be addressed 
in future publications. 

Note that our model is necessarily an oversimplified 
one. Firstly, it assumes that the state of passivation of 
the edge is solely determined by the concentration of Hi 



molecules present in the atmosphere of the experiment. 
This is not an entirely realistic assumption, since it is 
to be expected that the binding-unbinding of H atoms 
to an hexagon or heptagon (pentagons have no dangling- 
bonds to which H atoms can bind to) is also influenced 
by temperature and pressure. Taking such observation 
into account in our model would imply the introduction 
of a chemical potential regulating the chemical equilib- 
rium between the passivating atoms attached to the edge 
and those in the atmosphere surrounding the ribbon. 
Moreover, in the most general case, one would have also 
to allow the state of passivation of an hexagon or hep- 
tagon to be a statistical variable, since the bare exchange 
parameters between neighboring sites should depend on 
their state of passivation. This would imply the intro- 
duction of a Potts-like model with a higher number of 
colors (corresponding to both hydrogen-passivated and 
non-passivated edge polygons). 

Finally, it is to be expected that in an experiment, 
the edge may be passivated by other atomic or molecular 
species present in the gaseous environment surrounding 
the ribbon (namely oxygen, nitrogen, water, etc.) and 
not just by hydrogen. In order to take into account the 
presence of competing species, one would need to consider 
a Potts-like model with a yet higher number of colors, 
together with additional exchange parameters associated 
with the interaction between different kinds of passiva- 
tion between neighboring polygons, each of which would 
need to be computed from ab-initio simulations. More- 
over, one would have to introduce a chemical potential 
for each species, regulating the chemical equilibrium be- 
tween the passivating atoms of that species attached to 
the edge and those in the atmosphere surrounding the 
ribbon. This would make the model increasingly difficult 
to study using a simple analytical approach as the one 
presented above. 



III. CONCLUSION 

In this work, we have treated the graphene nano-ribbon 
edge as a one-dimensional system and introduced a 3- 
color Potts-like model for the zigzag edge reconstructed 
with Stone-Wales defects. We have shown how to extract 
the effective parameters that describe the model from ab- 
initio calculations and how to use these numerical val- 
ues to determine the temperature dependence of the de- 
fect size and concentration at the edge. We concluded 
that according to the simple model presented, both the 
totally passivated and the totally non-passivated zigzag 
edges always present a finite concentration of defects at 
any finite temperature, as is to be expected for a one- 
dimensional system with short-range interactions. The 
room-temperature defects concentration is found to be 
exponentially dependent on the exchange parameters of 
the model. In addition, the density of defective domains 
and their size sharply decreases with decreasing temper- 
ature. We have also computed the domain distribution 
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size of domains of polarized and unpolarized spins and 
have concluded that these distributions do not have fat 
tails. 
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Appendix A: Exchange parameters values 

In order to obtain indicative values for the exchange 
parameters of the 3-color Potts-like model introduced 
above, we have not only used ab-initio results on 
non-passivated zigzag edges already published in the 
literature,^ but we have also performed ab-initio cal- 
culations on hydrogen-passivated edges. In Fig. [9] we 
plot the edge energies (relative to the energy of the pris- 
tine zigzag edge), obtained from ab-initio calculations of 
edge reconstructed zigzag ribbons with both hydrogen- 
passivated edges and non-passivated edges. 

From the ab-initio results summarized in Fig. |9| an d 
after employing the method presented in Section |II A| in 
order to compute the exchange parameters in both cases, 
we obtain the following values for the exchange parame- 
ters: 

• Hydrogen-passivated edge (case C\): We have 
performed ab-initio calculations for hydrogen- 
passivated ribbons with various defect concentra- 
tions (see Fig. [9]). 7 We have assumed the same 
value <r n = 0.01 eV for all uncertainties. 7 The val- 
ues obtained for the parameters are 7^ = (0.53 ± 
0.03) eV and 5 Cl = (0.65 ± 0.04) eV. 

• Non-passivated edge (case C 2 ) : The ab-initio results 
used in this case were extracted from the work of 
Huang et al.Wl We have assumed the same value 
a n = 0.01 eV for all uncertainties. 7 The values ob- 
tained for the parameters are 7 Ca = (0.03±0.02) eV 
and 6 C2 = (-0.49 ± 0.03) eV. 

In Fig. [9] we plot the results obtained from ab-initio 
calculations (isolated dots), compared with the polyno- 
mial interpolation of these results (dashed-lines) , which 
was obtained from a least-squares method, as described 
in section |II A| The fact that these curves are in good 
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Figure 9. Energies of the partially reconstructed edges, mea- 
sured relative to the pristine zz edge, as a function of the 
defect concentration. Dots represent the edge energies ob- 
tained from ab-initio calculations, while dashed lines corre- 
spond to a polynomial interpolation of the results obtained 
from the Potts-like model using the exchange parameters com- 
puted with the least squares method. The red squares and 
red dashed line represent the results obtained for hydrogen- 
passivated zigzag ribbons, whereas the blue circles and blue 
dashed line represent the edge energies of non-passivated 
zigzag ribbons. 



agreement with the ab-initio results justifies a posteriori 
the use of a Potts-like model with only nearest-neighbor 
interactions. 

Using the above calculated exchange parameters, the 
density of polarized spins, n po i, defined in Eq. d9j), the 

density of links between polarized spins, n_| , defined in 

Eq. (10), the density of polarized domains, rid, defined 



in Eq. (Ill and the average domain size of the minor 
domains, La v , defined in Eq. (21 1, acquire the form pre- 
sented in Fig. 10 

Using the exchange parameters calculated above, one 
concludes that the density of polarized domains in the 
non-passivated case at room temperature has a value 
(rid) — 6.65 x 10 -18 defects per unit cell of pristine zigzag 
edge (or (rid) — 5.41 x 10 -8 defects per meter). 7 In ad- 
dition, at room temperature, the unpolarized domains 
average domain size of Lq IF ~ 1 unit cells (or 
1.23 x 10~ 10 m). These small unpolarized do- 
mains are on average 1.85 x 10 7 m apart from each other. 

In the case of a hydrogen-passivated edge, the above 
calculated parameters give, at room temperature, a den- 
sity of polarized domains of (rid) — 1.50 x 10 -22 per unit 
cell of pristine zigzag edge (or (nd) — 1.22 x 10~ 14 defects 
per meter). The polarized domains have a mean size of 
IMF _ 1 unit cellg ( or imf _ i.23x 10- 10 m). The small 
domains of polarized spins are on average 8.20 x 10 13 m 
apart from each other. 

These results show that the ground states in both the 
non-passivated case (totally reconstructed edge) and the 
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written, using the cyclic invariance of the trace, as 



2N 



Tr(aT 2Ar ) 



Z 



2N 



(Bl) 



where T a p are the individual matrix elements of the 
transfer matrix in Eq. while Z 2 n = Tr(T 2Ar ) = 

X 2 ^ + Xq N + A 2 ^ is the partition function of the model 
and <t is the 3x3 matrix 



(B2) 




Figure 10. Plot of the four thermodynamic quantities intro- 
duced in the main text, for the values of the exchange pa- 
rameters obtained from the ab-initio results: 5c 1 = 0.66 and 
7d = 0.52 for the hydrogen-passivated edge (red curves); 
Sc 2 = —0.49 and 7c 2 = 0.03 for the non-passivated edge (blue 
curves). Panel (a) shows the density of polarized domains. 
Panel (6) shows the density of links between polarized spins. 
Panel (c) shows the density of polarized spins. Panel (d) shows 
the average minor domain size. The dashed green flat line 
represents the infinite temperature limit for both cases. 



hydrogen-passivated case (pristine zigzag edge) are very 
stable with respect to the effect of thermal disorder. 

However, we should note that in our model, the sensi- 
tivity of the thermodynamic quantities on the values of 
the exchange parameters is large. In order to illustrate 
such fact, consider for instance that the exchange param- 
eters were reduced to 1/4 of the values which we deter- 
mined above. This would imply that the room tempera- 
ture concentration of defects would be increased by sev- 
eral orders of magnitude: in the non-passivated case, it 
would be increased to (rid) — 7.01 x 10 5 defects per meter 
(the average distance between neighboring defects would 
be 1.43/im); in the hydrogen-passivated case, the defect 
concentration would be increased to (rid) — 4.82 x 10 4 
defects per meter (the average distance between neigh- 
boring defects would be 20.8 fj,m apart). Consequently, 
if more detailed ab-initio calculations were to give sig- 
nificantly smaller exchange parameters, this would imply 
that the edge ground state would be much less robust to 
the effect of thermal disorder. 



Eq. (Bl) shows that the magnetization is space- 



independent. Since the trace in Eq. ( Bl I is independent 



of the basis used for its calculation, we choose the one 
that diagonalizes T, 



|Ao> 



1 

7^ 



( i 
o 

V 



(B3) 



-1 



which is the eigenvector corresponding to the eigenvalue 
A = -er^ and 



|A±> 




(B4) 



which are the eigenvectors corresponding to the eigenval- 



ues A± = | 1 + e-P & ± V(l-e-^) 2 +8e- 2 ' 3 '/ 



a± is given by 



a± = 



A+-1 



v / 2[(A±-l) 2 +2e- 2 ^)] 1 /2 



where 



(B5) 



and where 0± = 1 — 2ot± (normalization condition). It 
can be easily checked that these three vectors form an 
orthonormal basis. Expressing the trace in terms of this 
basis, one obtains for the magnetization Eq. (Bll, the 
result 



1 

Z2N 



£ Af (\,\a\\,)=Q, 



(B6) 



/x=0,±l 



Appendix B: Correlation functions of the Potts-like 
model 

The computation of correlation functions of a one- 
dimensional Potts-model in a periodic system involves 
the computation of the trace of a string of operators^!. 
For instance, the magnetization of the system can be 



since (A M | a |A^) = for each one of the eigenvectors of 
T. This equality merely reflects the symmetry of the 
model with respect to an interchange of + with — spins 
that is present by construction. In order to infer the ex- 
istence of a phase transition at T — in the absence of 
a (infinitesimal) field that explicitly breaks this symme- 
try, one needs to consider the behavior of higher-order 
correlation functions. 
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The spin-spin correlation function 

(<Ti<T l+j ) = — 
^2N 
1 



' Oi<Ji+j ) is given by where a 2 is the matrix 



Tr(T 



2N- 



^2N 



E A f~ 



J jT J a) 

j K | (A p | a \X U ) | 2 , (B7) 



where we have used a representation of the unit-operator 
in terms of the eigenstates of T, on going from the first 
to the second line of Eq. (B7). At T ^ and in the 
thermodynamic limit N — > oo, the only term in the nu- 
merator of Eq. (B7) that survives, is the one with v = 0, 
fi = +1 and Z2N ~ X 2 ^ . 
since (A + | a |A ) = V2a+, 



Thus, we obtain in this case. 



2a\ 



Ao 



(B8) 



If T 7^ 0, ( (Ti(Ji + j ) — ^ if j — ^ 00, showing that the mag- 
netization of the system is zero at any finite temperature, 
as is to be expected for any system with Z2 symmetry in 
Id. At T — 0, one has to distinguish three cases: 6 > 0, 
in which case A + — > 1 and both Aq and A_ go to zero. In 



that case, Eq. ( B7 1 still holds and the ground state is sim- 
ply the 0000 . . . state, with no associated magnetization. 
If, on the other hand 5 < 0, A + — > 00, Xq — > —00 and 
A_ — s- 0. In that case, one has to consider again Eq. (B6 1, 
since the terms (A + |cr|Ao) and (Ao|<r|A + ) contribute 
equally to it. Thus, we obtain (<7^<Tj-|_j ) = (— 1) J , which 

shows that the anti-ferromagnetic states '. . . H 1 — . . .' 

and '. . . — I h • • • are the two degenerate ground-states. 

In this case, the system shows a transition to a finite 
(staggered) magnetization at zero temperature. Finally, 
if 6 = 0, X± —> 1, Ao — > —1 and all terms (A±|<t|Ao) 
and (Ao|<t|A±) contribute to Eq. (B7). We obtain 
( <Ti<Ti + j ) = I (— l)-' , which shows that there are three de- 
generate ground states '. . . 0000 '. . . H 1 — and 

— I h...'. One can also show that a phase transition 

is present when 6 < 0, if one writes equation Eq. ( B8 1 as 
( (7^+1 ) = 2a\ {-iy e-J'/C, where £ = 1/ \n{X + /JX |) 
is the correlation length of the model. If 8 > 0, £ = 
at T — and no phase transition occurs, but if 6 < 0, 
£ — > 00 at T = indicating the presence of a phase tran- 
sition. Note that the presence of a phase transition has 
at most a marginal effect on the results presented in the 
main text, since such a phase transition is due to the 
existence of a Z2 symmetry in the model, whereas the 
formation of minority domains of either 0's or H — relics 
on states not related by such a symmetry. 

One can also use the transfer matrix formalism to 
compute the probability ( <5 CTi(Ti+ ,_i ) that the spins at 
sites i and i 
5, 



j are anti-parallel. 



,-1 



1 



One uses the identity 
, which can be easily 



checked by substituting <Tj and cn+j by their values 0, ±1. 
Since we have already computed the spin-spin correla- 
tion function above, we are left with the computation of 
(£r 2 <7 2 +3 ). Following the same steps as above, we obtain 
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Z2N 



Y.K N ~ j K\{x^ 2 \x v ) 



(B9) 




(B10) 



At T ^ 0. 
ered in Eq. 
(A+I^IAo). 
pression Eq 

(5, 



the only terms that need to be consid- 



(B9l are those involving (A + |<r 2 |A + ) and 



Taking this into account, as well as the ex- 
one finally obtains for 



■><Ti<J i+j , 



(B8l for 



le result 



>i+3 . 



2ai + 2a 2 + a 2 _ 

A4 



An 

A, 



(Bll) 



Using Eq. (B8l with j 



in the expression for n c i — { o t i 

we obtain = 2a 2 A\ 
result given in Eq. ( 14 1 



and Eq. (Bll) with j 
2\ _ / x 



1 



section 
exactly 



II B 
the" 



) given m 
/A + ), which is 



Appendix C: The Domain Size Distribution 

In what follows, we will compute the domain size dis- 
tribution (DSD) of the 0-spins domains and of ±-spin 
domains. Throughout the computation, we will assume 
PBCs for the system. We will illustrate the computation 
of the DSD for the domains of 0-spins, pinpointing the 
differences with the computation of the DSD of domains 
of ±-spins. 

In the context of the exact calculation of the DSD, the 
thermal average of the distribution of the sizes of domains 
of 0-spins is defined by 



L 



2N 

E 

L=l 



LP (L), 



(CI) 



where Pq{L) is the domain size distribution of domains 
of 0-spins (thermal average of the fraction of domains of 
size L). This quantity is given by 



Po(L) = 



N d0 



(C2) 



where N% stands for the number of domains of 0-spins 
with size L, while N^o stands for the total number of 
0-spins domains regardless of their size. An analogous 
expression can be written for the DSD of ±-spins. 

We start by defining the operator that verifies in every 
possible way if the spin i is in a domain of spins with 
size L, 



L-l 



fi,L = E 



fe=0 



L-l 



°?-fe-l( ]^[ (l-of_ fc+7 ; y>, h+L 

7=0 



(C3) 
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The above definition is valid for cases where the do- 
main has a size L < 2N — 2, where 2N stands for the 
total number of spins in the one-dimensional chain. In 
the case where L — 2N — 1 and L — 2N, this definition 
is modified. It then reads 



Ji,L= 



i,L=2N-l 



fi 



L=2N 



2N-2 

£ 

fc=0 L 
2N-1 

n a 

7=0 



2N-2 

. n 

7=0 



-fc+7 



,(C4) 



(C5) 



We can analogously define an operator verifying in ev- 
ery possible way if the spin i is in a domain of spins ± 
with size L, namely f^ L . To do this, it suffices to substi- 
tute, in Eqs. (C3)-(C5l, the operators a 2 by 1 — a 2 . For 



L < 2N - 2, /fTreads 



L-l 



&-£ 



k=0 L 



(1 



i—k — 



L-l 

on* 

7=0 



2 

i — fc+7 



(1 - Of-fc+i) 



whereas, for L = 27V — 1 and L = 27V, we have 



fi,L=2N-l 



f^L=2N 



2N-2 

£ 

fc=0 
2JV-1 

n< 

7=0 



'i+7- 



-fc-1. 



2JV-2 
7=0 



fc+7 



(C6) 



>(C7) 



(C8) 



Given this, if for a given configuration of the edge, 
we want to count the number of unpolarized spins in 
domains of size L, N®, we just need to perform the sum of 
the operators ff L over every site in the one-dimensional 
chain. Explicitly, it reads 



2N 



Nl 



E/u 

i=l 



(C9) 



From this, we can easily extract the number of 0-spin 
domains with size L of a particular configuration. We 
have thus that the number of domains of 0-spins with 
size L is given by 



dL - L 



2N 

E 



f 

Jj,L 



(CIO) 



In addition, if we want to count all the domains of 0- 
spin, irrespective of their size, we just have to sum 
over all possible sizes L, 



2!\ 



N, 



./(I 



E< 



L=l 



2N 2N 

EE'" 

L=l i=l 



f 
~L 



(Gil) 



Instead of defining the operator counting the tota l 
number of domains of 0-spins as was done in Eq. (Cll I, 



we can use an equivalent and simpler expression for such 



an operator. It reads Ndo = N±o/2 + 6l ,2N, where N±q 
stands for the number of links between polarized and un- 
polarized spins, while the term 5l 2N accounts for the 
situation in which all the spins in the chain are unpolar- 
ized, in which case there are no links between polarized 
and unpolarized spins, but there is one domain of 0-spins 
occupying the entire chain. This operator can be written 
explicitly as 



2N 



N, 



dO 



2 i 2 ,^2 2 



2N 



n(!-^ 



(C12) 



Such a definition must be introduced, because the op- 
erator N±q/2 is not equivalent to the operator counting 
the number of domains in a one-dimensional spin chain. 
This operator, in fact, counts the number of links be- 
tween polarized and unpolarized spins (+0, —0, 0+ and 
0—) divided by two. Whenever the spin configuration is 
such that there are links between polarized and unpo- 
larized spins, this operator is equivalent to the operator 
giving the number of domains. However, when there are 
no links between polarized and unpolarized spins, this 
operator always yields 0, not being able to distinguish 
between the cases where all the spins are polarized (and 
thus the number of unpolarized spin domains is N^o = 0) 
and the case where all the spins are unpolarized (and thus 
the number of unpolarized spin domains is N^q = 1)- In 
order to account for these cases, the term Yli=i (l — °T) 
is added to the definition of N^o, giving 1 when the whole 
spin chain is unpolarized. 

We can write analogous equations to Eqs. ( C9 1-( C12 1 
for the case of polarized spins. The operator counting 
the number of polarized spins in domains of ±-spins, JVJ, 
reads 



2N 
i=l 



(C13) 



while the operator counting the number of domains (of 
±-spins) with size L, reads 



N ± = £tL_ = y t±L. 

dL L ^ L 

i=l 



(C14) 



The total number of ±-spins domains, irrespectively of 
their size, Nd±, reads 



2A" 



2N 2N f± 
Ji.L 



(C15) 



L=l L = l i=l 



which, analogously with what was done for Ndo, can be 
rewritten, reading 



^ 2JV 2JV 
N d ± = 2 E + ^+1 ~ °i a i+l + ^^i+l] + II a 



(C16) 
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Note that the term Y\f=i a 1 > evaluates to 1 when all the 
spins are polarized (forming a polarized spin domain with 
a length L = 2N) and to in all other cases. 

Given this, we can obtain the DSD of domains of 
0-spins, computing the thermal average of the ratio 
N% L /Ndo [see Eq. ( C2 1] . The computation of thermal 
averages of ratios can be performed using the following 
mathematical trick 



N, 



,/(> 



Nl L e~ uNd0 Au 



o 

oo 



(N^ e -uN io y Uj (C17) 



where in the last equality we assumed that the combined 
sum is convergent, regardless of the number of spins in 
the system. 

In the above thermal average, given by the sum over 
all configurations, we need to exclude the two config- 
urations with all spins polarized, since N® = and 
Ndo = in such case, yielding indeterminate terms to 
the sum. This is equivalent the computation of the con- 
ditioned probability of having a domain of unpolarized 
spins with a particular size L, given that there are do- 
mains of 0-spins in the one-dimensional chain. Excluding 
these terms changes the partition function, Z 2 Ni from 
Z 2 n = A+^ + X^ N + X 2 ^ to Z 2N = Z 2 n — 2e~ 2N P s — 
X 2 ^ — X 2N + A 2 ^. In addition, note that the sums over 
all the configurations must also exclude the terms asso- 
ciated with this configuration. 



We can thus rewrite Eq. ( C17 1 as 



NS 



N, 



40 



N e -l3E({a})~uN d0 



{a}' 



-du, (C18) 



J 2N 



where {a}' indicates that the sum is performed over all 
the configurations except the two configurati ons w ith all 
spins polarized. Note, however that in Eq. ( C18 1 , sum- 
ming over {a}' or over all the configurations, {a\, yields 
the same result, because the two configurations with all 
spins polarized contribute with AT9 = to the sum. 

The version of Eq. ( C18 I for domai ns of polarized spins 
is obtained by substituting in Eq. (C18) N® and Ndo 



by respectively, N d and Nd±, while Z' 2N = Z 2 n — 1 



X\ N 



X_ — 1, since in this case the configuration 
yielding N^± = is that with all the spins unpolarized. 
Explicitly, it reads 



-j3E({v})-uN di 



Z: 



-du, (C19) 



2N 



a. The exact expression of the DSD 

In order to obtain the exact expressio n of th e DSD, we 
need to compute the integrand of Eq. (C18). The sum 
over all configurations in Eq. ( C18 ) can still be performed 



using the transfer matrix formalism (see Appendix |Bj . 
However, we must use a modified transfer matrix, which 
absorbs the exponential of the number of domains N^o, 
appearing in Eq. ( |C18 1 , in the definition of the T-matrix 
given in Eq. (151. It reads 



T 



( 



\ 





-06 




(C20) 



where we have rescaled the exchange parameter 7 to 
7 = 7 + u/(2/3). Both the eigenvalues, A + , Ao, A_ and 
the eigenvectors, | Ao) , | Ao) , | Ao), of this T-matrix, have 
exactly the same form of those obtained for the T-matrix, 
with 7 substituted by 7. Now, both the eigenvalues and 
the eigenvectors depend on the integration variable u, 
through the rescaled exchange parameter 7. 

As mentioned above, this new T-matrix originates from 
the definition of N^o [see Eq. ( C12 1] . Note however, that 
there is a subtlety in the definition of the T-matrix in Eq. 



(C2CM. In fact, this transfer matrix absorbs not e 



-uN d 



-uS L 



but instead e ~ uN±o/2 into itself. The factor 
that also enters the definition of Ndo [see Eq. ( C12 1] , is 
not absorbed into T-matrix, because this term involves 
all the spins of the chain, which cannot be properly rep- 
resented using a nearest neighbor transfer matrix formal- 
ism. As a consequence, we must keep in mind that the 
results of the sum over all configurations in Eq 



(C18I 



m the 
cases where To = 2N. 

For domains of polarized spins, the sum over all config- 
urations in Eq. (C19 1, is still computed using a modified 
transfer matrix 
Eq 



This T-matrix is the same as that of 
(C20), defined for the case of unpolarized domains. 



However, if we remember the definition of the operator 
counting the number of polarized spins domains [see Eq. 



( C16 1] , Nd± = N±o/2 + 5z, ± 2jVj we readily conclude that 



now, our results will need to include an additional fac- 
tor of e~ u in the cases when L± = 2N and not when 
L = 2N. 

Given this, we can rewrite the integrand in Eq. (C18l 
for the DSD of unpolarized spins domains, using the 
transfer matrix formalism as, 



h{L) 



1 Z 2 n 



L 



27V 



7' 

n 2N 



-u5r 



(C21) 



where Z 2 n is the partition function associated with the 
T-matrix. In what concerns the computation of the DSD 
of domains of ±-spin domains, note that the integrand 



i n Eq . ( C19 1 , I±(L), is of the same form as Iq(L) on Eq. 



(|C2l|, but with f° L substituted by ff h , 



I±(L) 



1 Z. 



2N 



2N 



Z 2N L ,_, 



-u8 L 



(C22) 



where we should recall that the ZL N in Eq. ( C22 1 is dif- 



ferent from that appearing in Eq. ( C21 1 . In addition 
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note that the exponential term in Eq. (C22l refers to 



the configuration where all the spins of the ID chain are 
polarized, while that at Eq. (C21 1, refers to the configu- 



ration where all the spins are unpolarized. 

If we now define O(L) = ■ ■ ■ C i+L-i where & 



1 - of, we can, using Eqs. (C3)-(C5l, write Eq. (C21) 



I {L < 2N — 2) = 2N^- (Q(L + 2)) f 



Z' 



2 A 



I (L = 2N - 1) = 2N 



2(9(L + l)) f + (Q(L)) f (C23a) 
Z2N 



Z' 



(6(2AT)}, 



J 2N 

(®{2N-l)) f 



I (L = 2N) = 



^»(G(2N)) f e-\ 

^2N 



(C23b) 
(C23c) 



In the case of the domains of polarized spins, in ac- 
cordance with Eqs. (C6l-(C8l, the integrands are ob- 



tained from Eqs. (C23), just by substituting 0(L) by 



I ± (L < 2N — 2) = 2N ^ (T(L + 2)) 



I±(L =2N-1) 



I±(L = 2N) 



7' 

n 2N 



2(r(L + l)) f + (T(L)) ¥ (C24a) 

Z 2 N 



2N- 



7' 
n 2N 



(T(2N)) f 



(T(2N-l)) f 
Z2N 



Z, 



{T{2N)) f e 



(C24b) 
(C24c) 



2N 



Computing the correlation functions (<d(L)}^, using 
the transfer matrix formalism involves the computation 
of the following trace 



(Q(L)) f = 



1 



^2N 



Tr 



T 



■2N-L 



(en 1 



(C25) 



which yields the result (6(L)) f = F(1) L - 1 F(2N -L + 
1)/(Z 2N F{0) L ), where F(r) = S+/L AH - a^+X r + . The 
5± and j3± are the entries of the eigenvectors of T, |A±) 



[see Eq. (JB5J)]. Noting that F(l) = F(0) and using Eq. 
(B5), we can finally write (Q(L))f as 



{0{L))t z 2N v ^M)W ^' (C26) 

where p = 2N - (L - 1). 

In the computation of the DSD of the polarized spin 
domains, the {T(L))f appearing in Eqs. (C24) can be 
analogously computed and one obtains 



(T(L)) f 



e -08(L-l) (A+ _ 1)A P _ ( A _ _ 1)A P 

Z 2N \J (e _/35 - l) 2 +~8e-2^ 



\2N 



->2N 



(C27) 



where, again, p = 2N — (L — 1). 

The integrands in Eqs. ( C23 1 can be rewritten as 



I (L <2N-2) = 27V -J_ (w_ - W+) , (C28a) 
Z 2N \ / 

I (L = 2N-l) = 2N-? r - (V_ - Y+) , (C28b) 

^2N K ' 



I (L = 2N) = — 



(C28c) 



2N 



where W± = A^ 2 (A ± - 1)(A T - 1) 2 /(A + - A_), while 

Y ± = X^CXi - 1)(A T - 1)/(A+ - A_). 

For the polarized spins domains case, the integrands 
111 Eqs. ( |C24| can be rewritten as 



I ± {L<2N -2) = 2N— (W--W+), (C29a) 
I ± {L = ZN-\)=2N-£-[y--y+), (C29b) 

Zj 2N 



\2N 

I ± (L = 2N) = 2 -^e- 

^2N 



(C29c) 



where W± = e -^ L -^ X P ± 2 (X ± - 1)(A± - e -^)7(A+ - 



e - W L-i) A P-i (Ad 



l)(Ad 



A_), while y± 

e -^)/(A + _-A_). ^ 

Performing the integral over u in the expressions for 
Iq(L) as given in Eqs.( |C28 1, leaves us with the following 
expressions for the DSD of domains of unpolarized spins 
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N, 



ili) 



L<2N-2 



L=2N-1 



N d0 



2N 1 

v7 Om+2 
Zj 2N Z 



1 



(m + l)(m + 2) 
2iV 1 



L=2N 



Z' 2N 8 



7' ' 
"2N 



jc(G m+1 + G™ +1 ) - Vc 2 +d(G m+1 - G™ +1 ) - 2c(2(c+ 1) 
| (G m+2 + G™+ 2 ) - 2" l+2 (l + (c + l) m + 2 ) 



m+l 



G±+G\ -4 l + (c+l 



(C30a) 
(C30b) 
(C30c) 



with m = p - 2, G±=c + 2± ^c 2 + d, c = e"' 35 - 1 
and d = 8e _,S7 . Recall that, as the above equations refer 
to the computation of the DSD of domains of 0-spins, in 



r 



we have that Z' 2N = Ai — A 



2.Y 



Eqs. (jC30 

The DSD of domains of polarized spins, is 



given from Eqs.(C29), by 



+ X 2N . 
analogously 




N, 



d± 



L<2N-2 



L=2N-1 



N d± 



L=2N 



2N e-^(^-i) 

Om+2 
1 

(m + l)(m + 2) 

27V e -^( 2W - 2 ) 
Z 2iV 8 



2Af 

^2^ 



gm+2 + grn+2 ~j _ yra+2 f ]_ + (1 _ c)" i+2 ) 



2 2m+2 g 



(G^+G 2 ,) -4(l + (1-c) 2 ) 



(C31a) 
(C31b) 
(C31c) 



with c = 1 - e^ 135 and £7 ± = 2 - c ± Vc 2 + d. As Eqs. 
( C31 1 refer to the DSD of domains of ±-spins, the Z' 2N 
appearing in Eqs. ( JC3l) reads = A 2 ^ + A 2Ar + A 2 ^ - 
1. 



6. T/ie thermodynamic limit of the DSD 

In the thermodynamic limit (27V —¥ oo), when the 
domain size is much smaller than the size of the sys- 
tem (L -C 27V), we can perform the limit m — > +oo in 
Eq. ( |C30a| . Note that as G+ > G_, c + 1 > 1 and 
G+ > 2(c+ 1), then we have that the DSD of domains of 
unpolarized spins, in the thermodynamic limit, is given 



by Eq. (16). Analogously, the DSD of domains of polar- 



ized spins, also in the thermodynamic limit, is given by 



Eq. (151 



In Fig. [TT] we have plotted the DSD of domains of 
unpolarized spins at a non-passivated edge, as well as 
the DSD of domains of polarized spins at a hydrogen- 
passivated one, for two sets of particular values of the 
exchange parameters. In these plots we have used, the 
values for the exchange parameters 7 and 5 computed in 



Appendix [A] i.e. 7 = 0.03 and S — —0.49 for the unpas- 
sivated edge and 7 = 0.52 and 6 — 0.66 for the hydrogen- 
passivated edge. We readily see that both DSDs strongly 
decay with increasing size of the domain. Note the de- 
pendence of such decay with the temperature: the larger 
the temperature is, the smaller the decay rate is. 

A distribution P{x) is said have fat tails if it displays a 
slower decrease than the normal distribution, (or, alter- 
natively, if it decreases with a power of x) when x — > 00. 
As a consequence, the moments of a fat-tailed distribu- 
tion diverge above a given order, characteristic of that 
distribution, and thus its characteristic function is not 
analytical at the origin. 

Let us now consider the question of whether the DSDs 
computed above display fat tails in the thermodynamic 
limit (let us represent the DSD generically as V(L)). Its 
characteristic function is given by the discrete Fourier 
transform 



+00 



(C32) 



71=1 



The characteristic functions of Vo(L) and of V±(L) are 
geometric scries, and hence easily computable, after 



1G 




which we obtain Eq. (18 I and Eq. (17) 



J le-20 
- le-40 
J le-60 



Figure 11. Plots of the DSD for two sets of values of the 
exchange parameters 7 and 8. In all four panels, we present 
the DSD for six different temperatures: in red, T = lOOOif ; 
in black, T = 750A"; in green, T = 500K; in blue, T = 300A"; 
in yellow, T = 250A"; in brown, T = 150K. In panels (a) 
and (6) we plot the DSD of domains of unpolarized spins, in 
an unpassivated edge (7 = 0.03 and 8 = —0.49). Panel (a) 
has linear scale axis, while panel (6) has a logarithmic scale 
y-axis. In panels (c) and (d) we plot the DSD of domains of 
polarized spins, in a hydrogen-passivated edge (7 = 0.52 and 
8 = 0.66). Panel (c) has linear scale axis, while panel (d) has 
a logarithmic scale y-axis. 



From Eq. ( 18 1 , or its geometric series form, we con- 
clude that all the derivatives of Vq(w) exist at w = 0, if 
A+ > 1. If the unpolarized spins are the minority spins 
in the chain, 5 < and we have A+ > 1 + y/2 > 1 for 
every j3 > 0, thus Vo{L) has no fat-tails. 



Similarly, from Eq. ( |17| , when the polarized spins are 
the minority spins in the spin chain, 5 > and we con- 
clude that all the derivatives of V±(w) exist at w = 0, 
if X+e^ s > 1. This is verified for every (3 > 0, since 
X+e p5 > 1 + V2 > 1. Again, we conclude that V±(L) 
has no fat-tails. 

The first moment of Vq{L) gives us the mean size of 
the domains of unpolarized spins in the thermodynamic 
limit, Lq. which reads 



(C33) 



w=0 



In the same way, we can compute the mean size of the 
domains of polarized spins (in the thermodynamic limit), 
L±. Both these two quantities are written, respectively, 
in Eqs. (20 1 and (19). These two quantities are plotted 
in Fig. 
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